index

Data preparation

Code

library(dplyr)
library(ggplot2)

morphometric_data <- utils::read.csv(here::here("data", "morphometric_data.csv"), sep = ";", header = T, dec = ".")

morpho_data <- morphometric_data %>%
  select(-c(variable_abbreviation, variable_unit)) %>%
  t() %>%
  as.data.frame() %>%
  janitor::row_to_names(row_number = 1) %>%
  `rownames<-`(NULL)%>%
  # delete for now (n=1)
  filter(species!= "Diaphus_sp")

# replace empty value by NA 
morpho_data[morpho_data ==""] <- NA

# Numeric variables
morpho_data[, 4:23] <- sapply(morpho_data[, 4:23], as.numeric)

Data summary

Code

morpho_data_summary <-morpho_data %>%
  group_by(species) %>%
  count(species)

#htmltools::tagList(DT::datatable(morpho_data_summary))

Missing data

The operculum width was the variable with the highest number of missing data (n = 30). This represents 4% of the data

2 species with only NAs for an entire trait:

• Eurpharynx pelecanoides : 7 traits
• Melanostigma atlanticum : 2 traits

data imputation

• mice algorithm: n imputation = 5, n iterations = 50
• plot comparison for Operculum Width

Code

#select numeric variables for imputation 
original_data <- morpho_data %>%
  select(1:23)

imputation <-
  mice::mice(
    original_data,
    m = 5,
    maxit = 50,
    printFlag = F
  )

imputed_data <- mice::complete(imputation)

comparison <- tibble::tibble(species = original_data$species,
                             standard_length = original_data$standard_length,
                             Original = original_data$operculum_width,
                             Imputed = imputed_data$operculum_width) %>%
  arrange(species, standard_length)

ggplot(comparison, aes(standard_length, Imputed, col = is.na(Original))) +
  geom_point(size = 2, alpha=0.8) +
  facet_wrap(~species, scales = "free") +
  labs(col = "Imputed?") +
  scale_color_manual(values = c("grey60", "firebrick2")) +
  theme_bw()

Species * traits

calctulate functional traits

Code

# calculate functional numeric traits
numeric_traits <- imputed_data %>%
  na.omit() %>%
  select(-individual_code) %>%
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume
  ) %>%
  group_by(species) %>%
  summarise(across(everything(), mean, na.rm = TRUE)) %>%
  arrange(species)

# categorical traits for species without NA
cat_morpho <- morpho_data %>%
  select(
    species,
    ventral_photophores,
    gland_head,
    chin_barbel,
    small_teeth,
    large_teeth,
    fang_teeth,
    retractable_teeth,
    internal_teeth,
    gill_raker_types,
    oral_gape_axis
  ) %>%
    na.omit() %>%
  distinct() %>%
  arrange(species)

# combined the two data frames
fish_traits <- numeric_traits %>%
  inner_join(cat_morpho, by = "species") %>%
  arrange(species) %>% 
  tibble::column_to_rownames("species")%>%
  # assign trait type 
  # as.factor for qualitative traits
  mutate_if(is.character, as.factor)%>%
  # as.ordered for ordered variables
  mutate_at(c("gill_raker_types", "oral_gape_axis"), as.ordered)

Traits density distribution

Code

density_plot_traits <- fish_traits[,1:16] %>% 
    tibble::rownames_to_column(var="species") %>% 
    tidyr::pivot_longer(!species, names_to = "traits", values_to = "values")

ggplot(density_plot_traits , aes(values)) +
  geom_histogram(bins = 10,
                 color = "darkgrey",
                 fill = "lightgray") +
  facet_wrap(~ traits, scales = "free") +
  theme_minimal()+
  theme(strip.text.x = element_text(size = 11, face = "bold"))

Code

## Display the table ----
htmltools::tagList(DT::datatable(fish_traits))

traits correlation

Code

M <-cor(numeric_traits[, c(-1)])

ggcorrplot::ggcorrplot(M, hc.order = TRUE, type = "lower",
                       lab = TRUE, tl.cex = 9, lab_size = 3)

Code

ggsave("corrplot.png", path = "figures")

Code

# list of species 
sp_names <- c(rownames(fish_traits), "Nannobrachium_atrum", "Cyclothone", "Stomias_boa_boa")

# taxonomic_families
taxonomic_families <- sp_names %>%
  as.data.frame() %>%
  `colnames<-`("species") %>% 
  mutate(
    family = case_when(
      species %in%
        c(
          "Benthosema_glaciale",
          "Ceratoscopelus_maderensis",
          "Diaphus_metopoclampus",
          "Lampanyctus_ater",
          "Lampanyctus_crocodilus",
          "Lampanyctus_macdonaldi",
          "Lobianchia_gemellarii",
          "Myctophum_punctatum",
          "Notoscopelus_bolini",
          "Notoscopelus_kroyeri",
          "Bolinichthys_supralateralis"
        ) ~ "Myctophidae",
      species %in% c(
        "Borostomias_antarcticus",
        "Chauliodus_sloani",
        "Malacosteus_niger",
        "Melanostomias_bartonbeani",
        "Stomias_boa"
      ) ~ "Stomiidae",
      species %in% c(
        "Holtbyrnia_anomala",
        "Holtbyrnia_macrops",
        "Maulisia_argipalla",
        "Maulisia_mauli",
        "Maulisia_microlepis",
        "Normichthys_operosus",
        "Searsia_koefoedi",
        "Sagamichthys_schnakenbecki"
      ) ~ "Platytroctidae",
      species %in% c("Sigmops_bathyphilus",
                     "Gonostoma_elongatum") ~ "Gonostomatidae",
      species %in% c(
        "Argyropelecus_hemigymnus",
        "Maurolicus_muelleri",
        "Argyropelecus_olfersii"
      ) ~ "Sternoptychidae",
      species == "Anoplogaster_cornuta" ~ "Anoplogastridae",
      species %in% c("Arctozenus_risso", "Paralepis_coregonoides") ~ "Paralepididae",
      species == "Bathylagus_euryops" ~ "Bathylagidae",
      species == "Cyclothone_sp" ~ "Gonostomatidae",
      species == "Derichthys_serpentinus" ~ "Derichthyidae",
      species == "Eurypharynx_pelecanoides" ~ "Eurypharyngidae",
      species == "Evermannella_balbo" ~ "Evermannellidae",
      species == "Lestidiops_sphyrenoides" ~ "Lestidiidae",
      species == "Melanostigma_atlanticum" ~ "Zoarcidae",
      species %in% c("Photostylus_pycnopterus",
                     "Xenodermichthys_copei") ~ "Alepocephalidae",
      species == "Serrivomer_beanii" ~ "Serrivomeridae"
    )
  )

Species * assemblages matrix

Number of trawl hauls per depth

• Epipelagic = 8
• Upper mesopelagic = 26
• Lower mesopelagic = 16
• Bathypelagic = 16

Code

# Metadata
metadata <-  utils::read.csv(here::here("data", "metadata.csv"), sep = ";", header = T, dec = ".")%>%
  # calculation of standardized biomass values (vertical  trawl opening * horizontal trawl opening * distance traveled)  
  mutate(volume_filtered = 24*58*distance)

ggplot(metadata, aes(x=depth))  +
  ylab ("Number of trawls")+
  xlab("Immersion depth (m)")+
  geom_histogram(binwidth=100, col="white", fill=alpha("black",0.55))+
  theme_light()+
  coord_flip()+ 
  scale_x_reverse()+
  labs(fill= "")+
  guides(fill="none")+
  scale_y_continuous(breaks = c(2,4,6,8,10,12))+
  theme(axis.text.x= element_text(size=12),
        axis.text.y= element_text(size=12),
        axis.title.y = element_text( size=12),
        axis.ticks = element_blank())

Code

# species biomass x depth  matrix 2002-2019 ----
# data_biomass_2002_2019 <- utils::read.csv(here::here("data", "data_evhoe_catch_2002_2019.csv"), sep = ";", header = T, dec = ".")%>%
#   replace(is.na(.), 0)%>%
#   as.data.frame()%>%
#   rename("species"="Code_Station")%>%
#   mutate(species= gsub(" ","_", species))%>%
#   filter(species%in%sp_names)%>%
#   t()%>%
#   as.data.frame()%>%
#   janitor::row_to_names(row_number = 1)%>%
#   mutate_if(is.character, as.numeric)%>%
#   tibble::rownames_to_column("Code_Station")%>%
#   filter(!Code_Station=="H0472")%>%
#   tidyr::pivot_longer(!Code_Station, names_to = "species", values_to = "Tot_V_HV")%>%
#   rename("Nom_Scientifique"="species")

data_biomass_2002_2019 <- utils::read.csv(here::here("data", "Biomass_density_deep_pelagic_fish _Bay_of_Biscay.csv"), 
                                          sep = ";", header = TRUE, dec = ",")%>%
  rename( Code_Station= "Event.label") %>%
  select(Code_Station, 9:290) %>%
  mutate_all(~replace(., is.na(.), 0)) %>%
  select(c(contains("..g."),Code_Station)) %>%
  select(Code_Station, everything()) %>%
  select(-contains("..g.m3.")) %>% 
  rename_all(~gsub("\\.\\.g\\.", "", .)) %>% 
  rename_all(~gsub("\\.", "_", .)) %>% 
  tidyr::pivot_longer(cols = -Code_Station, names_to = "Nom_Scientifique", values_to = "Tot_V_HV") %>% 
  filter(Nom_Scientifique%in%sp_names)%>%
  filter(!Code_Station=="H0472")%>%
  replace(is.na(.), 0)

# species biomass x depth  matrix 2021 ----
data_biomass_2021 <- utils::read.csv(here::here("data", "data_evhoe_catch_2021.csv"), sep = ";", header = T, dec = ".")%>%
  select(Nom_Scientifique, Tot_V_HV, Code_Station)%>%
  distinct()%>%
  mutate(Nom_Scientifique= gsub(" ","_",Nom_Scientifique))%>%
  filter(Nom_Scientifique%in%sp_names)

# species biomass x depth  matrix 2022 ----
data_biomass_2022 <- utils::read.csv(here::here("data", "data_evhoe_catch_2022.csv"), sep = ";", header = T, dec = ".")%>%
  select(Nom_Scientifique, Tot_V_HV, Code_Station)%>%
  distinct()%>%
  mutate(Nom_Scientifique= gsub(" ","_",Nom_Scientifique))%>%
  filter(Nom_Scientifique%in%sp_names)

# #merge all matrix ----
depth_fish_biomass <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
   rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered)%>%
  # divise biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
  select(species, depth, biomass_cpu)%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  select(-depth)%>%
  group_by(species, depth_layer)%>%
  mutate(biomass=sum(biomass_cpu))%>%
  select(-c(biomass_cpu))%>%
  distinct()%>%
  tidyr::pivot_wider(names_from = species, values_from = biomass)%>%
  replace(is.na(.), 0)%>%
  tibble::column_to_rownames(var = "depth_layer")%>%
  #change species name
  rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
  rename("Cyclothone_sp"="Cyclothone")%>%
  rename("Stomias_boa"="Stomias_boa_boa") %>%
  as.matrix()

• assemblages = depth layers
• biomass data = all EVHOE data 2002-2022 (in m^3)

Code

htmltools::tagList(DT::datatable(depth_fish_biomass))

Species depth distribution

Code

depth_distribution <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%  
  rename("station"="Code_Station")%>%
  left_join(metadata)%>%
  rename("species"="Nom_Scientifique")%>% 
  select(-c(station))%>%
  distinct()%>%
  filter(Tot_V_HV>0)%>%
  group_by(species)%>%
  mutate(biomass_tot = sum(Tot_V_HV))%>%
  ungroup()%>%
  group_by(species, depth)%>%
  mutate(biomass_depth = sum(Tot_V_HV))%>%
  select(species, biomass_depth, biomass_tot, depth)%>%
  distinct()%>%
  group_by(species, depth)%>%
  mutate(biomass_rel=biomass_depth/biomass_tot*100)%>%
  select(species, depth, biomass_rel)%>%
  mutate(biomass = as.integer(biomass_rel))%>%
  select(-biomass_rel)%>%
  tidyr::uncount(biomass)

# Order in function of median depth
depth_distribution$species = with(depth_distribution, reorder(species, depth, median))  

ggplot(depth_distribution,
       aes(x = depth, y = species, group = species)) + 
  ggridges::stat_density_ridges(geom="density_ridges", scale=1.5, alpha=0.6, rel_min_height = 0.005,
                                quantile_lines = TRUE, quantiles = 2, size = 0.4, col= "gray30")+
  theme_bw()+
  ylab(label = "")+ xlab("Immersion depth (m)")+
  theme(axis.text.y = element_text(size=13),
        axis.text.x = element_text(face="italic", size=10, angle=80,vjust = 0.5, hjust=0),
        axis.title.x = element_text(size=13),
        axis.title.y = element_text(size=13))+
  xlim(0, 2000)+
  scale_y_discrete(position = "right")+
  scale_x_reverse()+
  coord_flip()+
  guides(fill="none", col="none" ,alpha="none")

Code

ggsave("depth_distribution.png", path = "figures", height = 8, width = 12)

Traits types

The first column contains traits name. The second column contains traits type following this code:

• N: nominal trait (factor variable)

• O: ordinal traits (ordered variable)

• Q: quantitative traits (numeric values)

• 1/3 of the traits are nominal (traits related to teeth and photophores), need to give them different weights so as not to overestimate functional diversity?

Code

fish_traits_cat <- utils::read.csv(here::here("data", "fish_traits_cat.csv"), sep = ";", header = T, dec = ".")
htmltools::tagList(DT::datatable(fish_traits_cat))

1. CWM

Community Weighted Mean : somme de l’abondance relative d’une espèce x valeur du trait

• trait quantittatif : valeur moyenne du trait si on prend un individu au hasard dans l’assemblage

• trait catégoriel : proportion des espèces possédant ce trait, une valeur élevée peut indiquer soit qu’un grand nombre possèdent se trait ou que l’espèce avec la plus forte abondance relative possède ce trait

• données centrées-réduites

Code

# spxtraits.matrix ----
spxcom.matrix <-  depth_fish_biomass %>% 
  t() %>% 
  as.data.frame() %>% 
  relocate("Epipelagic", "Upper mesopelagic", "Lower mesopelagic","Bathypelagic" ) %>% 
  tibble::rownames_to_column("species") %>% 
  arrange(species) %>% 
  tibble::column_to_rownames("species") %>% 
  as.matrix()

library(dplyr)

spxtraits.matrix <- fish_traits %>%
  mutate(across(16:23, ~ case_when(. == "P" ~ 1, 
                                   . == "A" ~ 0, 
                                   TRUE ~ as.numeric(.))),
         across(24, ~ case_when(. == "A" ~ 1, 
                                . == "B" ~ 2, 
                                . == "C" ~ 3, 
                                TRUE ~ as.numeric(.)))) %>% 
  select(-c(gill_raker_types, oral_gape_axis)) %>% 
  as.matrix()

# Remove the "[,1]" suffix from column names
names(spxtraits.matrix) <- gsub("[,1]", "", names(spxtraits.matrix))

#check rownames
#rownames(spxtraits.matrix) == rownames(spxcom.matrix)

result_CWM <- FD::functcomp(spxtraits.matrix, t(spxcom.matrix)) 
#FD::functcomp(spxtraits.matrix, t(spxcom.matrix), CWM.type = "all")

#  Calculate Total biomass
total_biomass <- colSums(spxcom.matrix)

#  Calculate Relative biomass
sp_rel_biomass <- t(spxcom.matrix) / total_biomass

# Transpose the Relative biomass Matrix for Display
t_sp_rel_biomass <- t(sp_rel_biomass)

total_sum <- colSums(t(sp_rel_biomass))

# Initialize an empty data frame to store results
CWM_df <- data.frame(
  depth_layer = character(),
  trait = character(),
  total_sum = numeric(),
  weighted_mean = numeric(),
  stringsAsFactors = FALSE
)

# Loop through each trait
for (trait in colnames(spxtraits.matrix)) {
  # Calculate the weighted sum for the current trait
  weighted_sum <- colSums(t_sp_rel_biomass * spxtraits.matrix[, trait])
  
  # Create a data frame for the current trait
  trait_df <- data.frame(
    depth_layer = colnames(t_sp_rel_biomass),
    trait = trait,
    total_sum = total_sum,
    weighted_mean = weighted_sum
  )
  
  # Append results to the main data frame
  CWM_df <- rbind(CWM_df, trait_df)
}

CWM_df <- CWM_df %>% 
  mutate(traits_names= gsub("_"," ", trait)) 

biomass <- t_sp_rel_biomass %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column(var = "species") %>% 
  tidyr::pivot_longer(!species, names_to = "depth_layer", values_to = "biomass") 

CWM_df <- CWM_df %>% 
  filter(!trait%in% c("chin_barbel",
                      "dorsal_fin_insertion", "eye_position",
                      "gland_head", "oral_gape_position", "oral_gape_shape",
                      "pectoral_fin_position", "retractable_teeth",
                      "transversal_shape")) %>% 
  mutate(traits_names= gsub("_"," ", trait)) 

CWM_df$depth_layer <- factor(CWM_df$depth_layer, 
                             levels = c("Epipelagic", "Upper mesopelagic",
                                        "Lower mesopelagic", "Bathypelagic"))

CWM_df$traits_names <- factor(CWM_df$traits_names, 
                              levels = c( "caudal throttle width", "oral gape surface",
                                         "large teeth", "eye size",
                                         "orbital length","small teeth",
                                         "internal teeth", "lower jaw length",
                                         "pectoral fin insertion", "fang teeth",
                                         "operculum volume", "ventral photophores",
                                         "gill outflow", "head length","body depth"))


ggplot(CWM_df, aes(x = depth_layer, y = weighted_mean, group = depth_layer, color = depth_layer)) +
  geom_point(position = position_dodge(width = 0.75), size = 3) +
  facet_wrap(~traits_names, scales = "free", ncol = 3) + 
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  labs(x = "",
       y = "Community Weighted Mean ") +
  guides(col="none")+
  theme_light()+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text(size =14, face="bold"),  
        legend.title = element_text(size =11),  
        legend.text = element_text(size =11), 
        axis.title.y = element_text(size=11),
        axis.text.y = element_text(size=12),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 15, face = "bold", color = "black"))

Code

ggsave("CWM.png", path = "figures", dpi = 700, height = 10, width = 10)

CWM bootstrap

• non parametric

Code

library(traitstrap)

# Trait 
trait_boot <- morpho_data%>% 
  inner_join(metadata) %>% 
  select(-c(individual_code, years, longitude_start,
            latitude_start, longitude_end, longitude_end,
            volume_filtered, distance)) %>% 
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    depth,
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume,
    ventral_photophores, 
    gland_head,
    chin_barbel, 
    small_teeth, 
    large_teeth, 
    fang_teeth, 
    retractable_teeth, 
    internal_teeth
  ) %>%
  mutate_at(vars(ventral_photophores, 
                 gland_head,
                 chin_barbel, 
                 small_teeth, 
                 large_teeth, 
                 fang_teeth, 
                 retractable_teeth, 
                 internal_teeth), 
            funs(ifelse(. == "P", 1, ifelse(. == "A", 0, .)))) %>% 
  mutate(across(all_of(c("ventral_photophores", 
                         "gland_head",
                         "chin_barbel", 
                         "small_teeth", 
                         "large_teeth", 
                         "fang_teeth", 
                         "retractable_teeth", 
                         "internal_teeth")), as.numeric)) %>% 
  tidyr::pivot_longer(!c(species,depth), names_to = "trait", values_to = "values")%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

# Community 
community <-  rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered)%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(Tot_V_HV))%>%
  select(-c(Tot_V_HV))%>%
  distinct()%>%
  select(-c(volume_filtered)) %>% 
  filter(biomass>0) %>% 
  mutate(species = case_when(
    species == "Nannobrachium_atrum"~"Lampanyctus_ater",
    species == "Cyclothone"~"Cyclothone_sp",
    species == "Stomias_boa_boa"~"Stomias_boa",
    TRUE ~ species
  )) 

trait_filling <- trait_fill(
  # input data (mandatory)
  comm = community,
  traits = trait_boot,
  
  # specifies columns in your data (mandatory)
  abundance_col = "biomass",
  taxon_col = "species",
  trait_col = "trait",
  value_col = "values",
  
  # specifies sampling hierarchy
  scale_hierarchy = c("depth_layer", "depth"),
  
  # min number of samples
  min_n_in_sample = 9
)

# run nonparametric bootstrapping
np_bootstrapped_moments <- trait_np_bootstrap(
  trait_filling, 
  nrep = 100
)

np_bootstrapped_moments_plot <- np_bootstrapped_moments %>% 
  mutate(trait= gsub("_"," ", trait))%>%
  filter(
    trait %in% c(
      "caudal throttle width",
      "large teeth",
      "oral gape surface",
      "gill outflow",
      "internal teeth",
      "oral gape shape",
      "small teeth",
      "orbital length",
      "transversal shape",
      "operculum volume",
      "body depth",
      "eye size"
    )
  )

np_bootstrapped_moments_plot$trait <- factor(
  np_bootstrapped_moments_plot$trait,
  levels = c(
    "caudal throttle width",
    "large teeth",
    "oral gape surface",
    "gill outflow",
    "internal teeth",
    "oral gape shape",
    "small teeth",
    "orbital length",
    "transversal shape",
    "operculum volume",
    "body depth",
    "eye size"
  )
)

np_bootstrapped_moments_plot$depth_layer <- factor(np_bootstrapped_moments_plot$depth_layer, 
                                              levels = c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic"))

# Mean
ggplot(np_bootstrapped_moments_plot, aes(x=depth_layer, y=mean)) +
  geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.1) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  labs(col="Depth layer", fill="Depth layer", y="Bootstrapped CWM")+
  facet_wrap(~trait, scales = "free")+
  theme_light()+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text(size =14, face="bold"),  
        legend.title = element_text(size =11),  
        legend.text = element_text(size =11), 
        axis.title.y = element_text(size=11),
        axis.text.y = element_text(size=12),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 10, color = "black"))

Code

ggsave("CWM_boot.png", path = "figures/additional_analyses", dpi = 700, height = 8, width = 12)

Summarizes bootstrapping output

Code

# summarizes bootstrapping output
sum_boot_moment <- trait_summarise_boot_moments(
  np_bootstrapped_moments
)

htmltools::tagList(DT::datatable(sum_boot_moment))

kruskal test beetween depth layers

Code

# Function to apply Kruskal-Wallis test and effect size calculation to each group
kruskal_with_effsize <- function(data, trait_name) {
  kruskal_res <- rstatix::kruskal_test(data, mean ~ depth_layer)
  effsize_res <- rstatix::kruskal_effsize(data, mean ~ depth_layer)
  combined_res <- cbind(trait = trait_name, kruskal_res, effsize_res)
  return(combined_res)
}

# Perform the Kruskal-Wallis test and calculate the eta² statistic for each trait
res.kruskal <- np_bootstrapped_moments_plot %>%
  group_by(trait) %>%
  group_map(~ kruskal_with_effsize(.x, .y$trait)) %>%
  bind_rows()

# View the results
htmltools::tagList(DT::datatable(res.kruskal))

Extracting raw distributions

Code

# run nonparametric bootstrapping
raw_dist_np <- trait_np_bootstrap(
  filled_traits = trait_filling,
  raw = TRUE
)

raw_dist_np$depth_layer <- factor(
  raw_dist_np$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
) 

raw_dist_np <- raw_dist_np%>% 
    mutate(trait= gsub("_"," ", trait)) %>% 
  filter(!trait%in% c("chin barbel", "fang teeth", 
                      "gland head", "internal teeth",
                      "large teeth", "retractable teeth", 
                      "small teeth", "ventral photophores"))

ggplot(raw_dist_np, aes(x = log(values), fill = depth_layer, col= depth_layer)) +
  geom_density(alpha = 0.2) +
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  labs(x = " log(Trait value)", col="Depth layer", fill="Depth layer") +
  facet_wrap(facets = vars(trait), scales = "free")+
  theme_light()+
  theme(strip.text.x = element_text(size = 10, face = "bold", color="black"),
    strip.background = element_rect(fill = "white"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10))

Code

ggsave("raw_traits_distribution.png", path = "figures", dpi = 700, height = 9, width = 12)

CWM - linear relationships

• selection R² < 0.6

Code

linear_relation <- np_bootstrapped_moments %>%
    mutate(trait= gsub("_"," ", trait))%>%
  filter(
    trait %in% c(
      "caudal throttle width",
      "large teeth",
      "oral gape surface",
      "gill outflow",
      "internal teeth",
      "oral gape shape",
      "small teeth",
      "orbital length",
      "transversal shape",
      "operculum volume",
      "body depth",
      "eye size"
    )
  )

linear_relation$trait <- factor(
  linear_relation$trait,
  levels = c(
    "caudal throttle width",
    "large teeth",
    "oral gape surface",
    "gill outflow",
    "internal teeth",
    "oral gape shape",
    "small teeth",
    "orbital length",
    "transversal shape",
    "operculum volume",
    "body depth",
    "eye size"
  )
)

ggplot(linear_relation, aes(x=depth, y=mean)) +
  geom_point(alpha=0.1, col="grey", size=0.8) +
  geom_smooth(method = "lm",col="#9565E5", se=T, linewidth=1, fill="#9565E5", alpha=0.5)+
  facet_wrap(~trait, scales = "free", ncol = 4)+
  theme_light()+
  ggpmisc::stat_poly_eq(formula = y ~ x, 
                        aes(label = paste(..rr.label.., ..p.value.label.. 
                                          , ..n.label..,sep = "*`,`~")),
                        parse = TRUE,
                        size=2.5,
                        label.x.npc = "right",
                        label.y.npc = "bottom",
                        vstep = -0.0005)+ 
  xlab("Trawling depth (m)")+
  ylab("Bootstrapped CWM")+
  theme(axis.title =  element_text(size =11), 
        strip.text = element_text(size =11, face="bold"), 
        axis.text= element_text(size=9),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 10, color = "black"))

Code

ggsave("CWM_boot_linear.png", path = "figures/additional_analyses", dpi = 700, height = 6, width = 10)

2. Build a functional space using the mFD package

2.1 Compute data summaries

Code

## Summary of the assemblages * species data.frame ----
asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)
asb_sp_fish_occ  <- asb_sp_fish_summ$"asb_sp_occ"

htmltools::tagList(DT::datatable(asb_sp_fish_occ))

2.2 Computing distances between species based on functional traits

• We have non-continuous traits so we use the Gower distance (metric = “gower”) as this method allows traits weighting.
• scale_euclid = TRUE

Code

sp_dist_fish <- mFD::funct.dist(
  sp_tr         = fish_traits,
  tr_cat        = fish_traits_cat,
  metric        = "gower",
  scale_euclid  = "scale_center",
  ordinal_var   = "classic",
  weight_type   = "equal",
  stop_if_NA    = TRUE)

## Output of the function mFD::funct.dist() ----
#round(sp_dist_fish, 3)

2.3 Building functional spaces and chosing the best one

2.3.1 Computing several multimensional functional spaces and assessing their quality

• mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).

Code

fspaces_quality_fish <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 10,
  deviation_weighting = "absolute",
  fdist_scaling       = FALSE,
  fdendro             = "average")

## Quality metrics of functional spaces ----
round(fspaces_quality_fish$"quality_fspaces", 3)

               mad
pcoa_1d      0.143
pcoa_2d      0.079
pcoa_3d      0.050
pcoa_4d      0.030
pcoa_5d      0.023
pcoa_6d      0.017
pcoa_7d      0.014
pcoa_8d      0.015
pcoa_9d      0.017
pcoa_10d     0.019
tree_average 0.041

The space with the best quality has the lowest quality metric. 5-D space good ?

Variance explained by each axis

Code

# Extract eigenvalues information
eigenvalues_info <- fspaces_quality_fish$"details_fspaces"$"pc_eigenvalues"

# Create a dataframe to store the results
variance_df <- data.frame(
  PC = c("PC1", "PC2", "PC3", "PC4"),
  VarianceExplained = c(
    eigenvalues_info[1, "Cum_corr_eig"] * 100,
    (eigenvalues_info[2, "Cum_corr_eig"] - eigenvalues_info[1, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[3, "Cum_corr_eig"] - eigenvalues_info[2, "Cum_corr_eig"]) * 100,
    (eigenvalues_info[4, "Cum_corr_eig"] - eigenvalues_info[3, "Cum_corr_eig"]) * 100
  )
)

htmltools::tagList(DT::datatable(variance_df))

2.3.2 Illustrating the quality of the functional spaces

Code

mFD::quality.fspaces.plot(
  fspaces_quality            = fspaces_quality_fish,
  quality_metric             = "mad",
  fspaces_plot               = c("tree_average", "pcoa_2d", "pcoa_3d", 
                                 "pcoa_4d", "pcoa_5d", "pcoa_6d"),
  name_file                  = NULL,
  range_dist                 = NULL,
  range_dev                  = NULL,
  range_qdev                 = NULL,
  gradient_deviation         = c(neg = "darkblue", nul = "grey80", pos = "darkred"),
  gradient_deviation_quality = c(low = "yellow", high = "red"),
  x_lab                      = "Trait-based distance")

This function generates a figure with three panels (in rows) for each selected functional space (in columns). Each column represents a functional space, the value of the quality metric is written on the top of each column. The x-axis of all panels represents trait-based distances. The y-axis is different for each row:

• on the first (top) row, the y-axis represents species functional distances in the multidimensional space. Thus, the closer species are to the 1:1 line, the better distances in the functional space fit trait-based ones.
• on the second row, the y-axis shows the raw deviation of species distances in the functional space compared to trait-based distances. Thus, the raw deviation reflects the distance to the horizontal line.
• on the third row (bottom), the y-axis shows the absolute or squared deviation of the (“scaled”) distance in the functional space. It is the deviation that is taken into account for computing the quality metric.

2.3.3 Testing the correlation between functional axes and traits

Code

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

# As we have 26 traits we have to split the df to see correlation between functional axes and traits 
# first set ----
fish_traits_1 <- fish_traits%>%
  select(1:9)

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_1, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value" < 0.05), ]

                trait axis         test stat value p.value
1            eye_size  PC1 Linear Model   r2 0.166  0.0081
2            eye_size  PC2 Linear Model   r2 0.160  0.0096
3            eye_size  PC3 Linear Model   r2 0.140  0.0158
4            eye_size  PC4 Linear Model   r2 0.139  0.0164
5      orbital_length  PC1 Linear Model   r2 0.398  0.0000
6      orbital_length  PC2 Linear Model   r2 0.114  0.0306
10       gill_outflow  PC2 Linear Model   r2 0.520  0.0000
13  oral_gape_surface  PC1 Linear Model   r2 0.136  0.0176
14  oral_gape_surface  PC2 Linear Model   r2 0.499  0.0000
18    oral_gape_shape  PC2 Linear Model   r2 0.141  0.0157
24 oral_gape_position  PC4 Linear Model   r2 0.200  0.0034
26   lower_jaw_length  PC2 Linear Model   r2 0.694  0.0000
29        head_length  PC1 Linear Model   r2 0.315  0.0001
30        head_length  PC2 Linear Model   r2 0.400  0.0000
34         body_depth  PC2 Linear Model   r2 0.360  0.0000
35         body_depth  PC3 Linear Model   r2 0.196  0.0038

Code

## Plot ----
fish_tr_faxes$"tr_faxes_plot"

Code

# second set ----
fish_traits_2 <- fish_traits%>%
  select(10:18)

fish_tr_faxes_2 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_2, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value" < 0.05), ]

                    trait axis           test stat value p.value
2   pectoral_fin_position  PC2   Linear Model   r2 0.268  0.0005
5  pectoral_fin_insertion  PC1   Linear Model   r2 0.272  0.0005
6  pectoral_fin_insertion  PC2   Linear Model   r2 0.412  0.0000
9       transversal_shape  PC1   Linear Model   r2 0.114  0.0309
11      transversal_shape  PC3   Linear Model   r2 0.228  0.0016
14  caudal_throttle_width  PC2   Linear Model   r2 0.402  0.0000
19   dorsal_fin_insertion  PC3   Linear Model   r2 0.170  0.0073
23           eye_position  PC3   Linear Model   r2 0.200  0.0034
26       operculum_volume  PC2   Linear Model   r2 0.126  0.0227
32    ventral_photophores  PC4 Kruskal-Wallis eta2 0.590  0.0000
33             gland_head  PC1 Kruskal-Wallis eta2 0.245  0.0011

Code

## Plot ----
fish_tr_faxes_2$"tr_faxes_plot"

Code

# third set ----
fish_traits_3 <- fish_traits%>%
  select(19:25)

fish_tr_faxes_3 <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits_3, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = T)

## Print traits with significant effect ----
fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value" < 0.05), ]

              trait axis           test stat value p.value
1       chin_barbel  PC1 Kruskal-Wallis eta2 0.142  0.0107
4       chin_barbel  PC4 Kruskal-Wallis eta2 0.183  0.0043
5       small_teeth  PC1 Kruskal-Wallis eta2 0.323  0.0002
9       large_teeth  PC1 Kruskal-Wallis eta2 0.467  0.0000
10      large_teeth  PC2 Kruskal-Wallis eta2 0.205  0.0027
13       fang_teeth  PC1 Kruskal-Wallis eta2 0.410  0.0000
21   internal_teeth  PC1 Kruskal-Wallis eta2 0.073  0.0499
23   internal_teeth  PC3 Kruskal-Wallis eta2 0.626  0.0000
25 gill_raker_types  PC1 Kruskal-Wallis eta2 0.327  0.0007
26 gill_raker_types  PC2 Kruskal-Wallis eta2 0.364  0.0004

Code

## Plot ----
fish_tr_faxes_3$"tr_faxes_plot"

Summary of traits with a significant effect

Code

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

fish_tr_faxes <- mFD::traits.faxes.cor(
  sp_tr          = fish_traits, 
  sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], 
  plot           = F)

## Print traits with significant effect ----
traits_effect <- fish_tr_faxes[which(fish_tr_faxes$p.value< 0.05),] %>% 
  as.data.frame() %>% 
  arrange(axis, desc(value))

htmltools::tagList(DT::datatable(traits_effect))

2.4 Plotting the selected functional space and position of species

Code

sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

big_plot <- mFD::funct.space.plot(
  sp_faxes_coord  = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
  faxes           = c("PC1", "PC2", "PC3", "PC4"),
  name_file       = NULL,
  faxes_nm        = NULL,
  range_faxes     = c(NA, NA),
  plot_ch         = TRUE,
  shape_pool      = 20,
  size_pool = 2.5,
  size_vert = 1.5,
  color_ch = "darkgrey",
  color_vert      = "black",
  fill_vert       = "black",
  color_pool = "grey",
  plot_vertices   = TRUE,
  check_input     = TRUE)

big_plot$"patchwork"

Code

ggsave("functional_space.png", path = "figures", dpi = 700, height = 9, width = 9)

Code

big_plot$PC1_PC2

Code

ggsave("functional_space_PC1_PC2.png", path = "figures", dpi = 700, height = 4, width = 4)

Code

big_plot$PC3_PC4

Code

ggsave("functional_space_PC3_PC4.png", path = "figures", dpi = 700, height = 4, width = 4)

3. Computing and plotting FD indices using the mFD package

3.1 Computing and plotting alpha FD indices

Code

alpha_fd_indices_fish <- mFD::alpha.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_w         = depth_fish_biomass,
  scaling          = TRUE,
  check_input      = TRUE,
  details_returned = TRUE)

The function has two main outputs:

• a data.frame gathering indices values in each assemblage (for FIde values, there are as many columns as there are axes to the studied functional space).

Code

fd_ind_values_fish <- alpha_fd_indices_fish$"functional_diversity_indices"
htmltools::tagList(DT::datatable(round(fd_ind_values_fish, 3)))

Code

fd_ind_values_fish_df <- as.data.frame(fd_ind_values_fish) %>% 
  tibble::rownames_to_column(var = "depth_layer") %>% 
  tidyr::pivot_longer(!depth_layer, names_to = "indices", values_to = "values" )

fd_ind_values_fish_df$depth_layer <- factor(fd_ind_values_fish_df$depth_layer, 
                                            levels = c("Epipelagic", "Upper mesopelagic",
                                                       "Lower mesopelagic", "Bathypelagic"))

fd_ind_values_fish_df$indices <- factor(fd_ind_values_fish_df$indices,
                                        levels = c("sp_richn", "fric", "fdis", "fdiv", 
                                        "feve", "fspe", "fide_PC1", "fide_PC2", "fide_PC3",
                                        "fide_PC4"),
                                        labels = c("Species richness", "Functional richness", "Functional dispersion", 
                                                   "Functional divergence", 
                                                   "Functional evenness", "Functional specialization",
                                                   "Functional identity PC1", "Functional identity PC2", 
                                                   "Functional identity PC3", "Functional identity PC4"))

fd_ind_values_fish_df2 <- fd_ind_values_fish_df %>% 
  filter(indices%in%c("Functional dispersion", "Functional richness",
                      "Functional divergence", "Functional specialization",
                      "Functional evenness"))

ggplot(fd_ind_values_fish_df2 , aes(x=depth_layer, y=values, fill=depth_layer)) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  geom_bar(stat="identity", position=position_dodge(), alpha=0.8)+
  facet_wrap(~indices, ncol=3)+
  theme_minimal()+
  labs(fill="Depth layer")+
  theme(axis.text.x = element_blank(),
        legend.text = element_text(size =10),
        axis.title.x = element_blank(), 
        strip.text = element_text(face="bold", size=11), 
        panel.border = element_blank())

Code

ggsave("indices.png", path = "figures", dpi = 700, width = 8)

Code

sp_richness <- fd_ind_values_fish_df %>% 
  filter(indices%in%c("Species richness"))

ggplot(sp_richness , aes(x=depth_layer, y=values, fill=depth_layer)) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  geom_bar(stat="identity", position=position_dodge(), alpha=0.8)+
  facet_wrap(~indices, ncol=3)+
  theme_minimal()+
 guides(fill="none")+
  theme(axis.text.x = element_blank(),
        legend.text = element_text(size =10),
        axis.title.x = element_blank(), 
        strip.text = element_text(face="bold", size=13), 
        panel.border = element_blank())

Code

ggsave("species_richness.png", path = "figures", dpi = 700, width = 3, height = 3)

functional identity

Code

fd_identity <- fd_ind_values_fish_df %>% 
  filter(indices%in%c("Functional identity PC1",
                      "Functional identity PC2", "Functional identity PC3",
                      "Functional identity PC4"))

ggplot(fd_identity, aes(x=depth_layer, y=values, fill=depth_layer)) +
  # paletteer::scale_fill_paletteer_d("ggsci::category20c_d3")+
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  geom_bar(stat="identity", position=position_dodge(), alpha=0.8)+
  facet_wrap(~indices, ncol=2)+
  theme_minimal()+
  guides(fill="none")+
  theme(axis.text.x = element_blank(),
        legend.text = element_text(size =10),
        axis.title.x = element_blank(), 
        strip.text = element_text(face="bold"))

Code

ggsave("fd_identity.png", path = "figures", dpi = 700, width = 6, height = 5)

• a details list of data.frames and lists gathering information such as coordinates of centroids, distances and identity of the nearest neighbour, distances to the centroid, etc. The user does not have to directly use it but it will be useful if FD indices are then plotted. It can be retrieved through:

Code

details_list_fish <- alpha_fd_indices_fish$"details"

Code

plots_alpha <- mFD::alpha.multidim.plot(
  output_alpha_fd_multidim = alpha_fd_indices_fish,
  plot_asb_nm              = c("Epipelagic", "Bathypelagic"),
  ind_nm                   = c("fdis", "fric", "fdiv", 
                              "fspe", "fide", "feve"),
  faxes                    = NULL,
  faxes_nm                 = NULL,
  range_faxes              = c(NA, NA),
  plot_sp_nm               = NULL,
  save_file                = FALSE,
  check_input              = TRUE) 

FRic Functional Richness

• the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.

• Only between epipelagic and bathypelagic layers

Code

plots_alpha$"fric"$"patchwork"

FDiv Functional Divergence

• the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage

Code

plots_alpha$"fdiv"$"patchwork"

FEve Functional Evenness

• the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.

Code

plots_alpha$"feve"$"patchwork"

FSpe Functional Specialization

• the biomass weighted mean distance to the mean position of species from the global pool (present in all assemblages).

Code

plots_alpha$"fspe"$"patchwork"

FDis Functional Dispersion

• the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.

Code

plots_alpha$"fdis"$"patchwork"

FIde Functional Identity

• the mean traits values for the assemblage. FIde is always computed when FDis is computed.

Code

plots_alpha$"fide"$"patchwork"

3.2.Computing and plotting beta FD indices

• The function returns a list containing:
• a dist object with beta indices values for each pair of assemblages:

Code

beta_fd_indices_fish <- mFD::beta.fd.multidim(
      sp_faxes_coord   = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")],
      asb_sp_occ       = asb_sp_fish_occ,
      check_input      = TRUE,
      beta_family      = c("Jaccard"),
      details_returned = TRUE)

Serial computing of convex hulls shaping assemblages with conv1

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
Serial computing of intersections between pairs of assemblages with inter_geom_coord

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |======================================================================| 100%

Code

head(beta_fd_indices_fish$"pairasb_fbd_indices", 10)

$jac_diss
                  Upper mesopelagic Bathypelagic Epipelagic
Bathypelagic              0.2597064                        
Epipelagic                0.5162017    0.6417795           
Lower mesopelagic         0.1632621    0.3583710  0.4417311

$jac_turn
                  Upper mesopelagic Bathypelagic   Epipelagic
Bathypelagic           0.000000e+00                          
Epipelagic             2.545169e-04 0.000000e+00             
Lower mesopelagic      3.697387e-02 6.636871e-16 6.778430e-05

$jac_nest
                  Upper mesopelagic Bathypelagic Epipelagic
Bathypelagic              0.2597064                        
Epipelagic                0.5159472    0.6417795           
Lower mesopelagic         0.1262882    0.3583710  0.4416633

• a vector containing the FRic value for each assemblage retrieved through the details_beta list:

Code

beta_fd_indices_fish$"details"$"asb_FRic"

Upper mesopelagic      Bathypelagic        Epipelagic Lower mesopelagic 
        0.7402936         1.0000000         0.3582205         0.6416290 

• a list of vectors containing names of species being vertices of the convex hull for each assemblage retrieved through the details_beta list:

Code

beta_fd_indices_fish$"details"$"asb_vertices"

$`Upper mesopelagic`
 [1] "Stomias_boa"                "Malacosteus_niger"         
 [3] "Melanostigma_atlanticum"    "Argyropelecus_olfersii"    
 [5] "Xenodermichthys_copei"      "Serrivomer_beanii"         
 [7] "Arctozenus_risso"           "Borostomias_antarcticus"   
 [9] "Maurolicus_muelleri"        "Sagamichthys_schnakenbecki"
[11] "Benthosema_glaciale"        "Lestidiops_sphyrenoides"   
[13] "Lobianchia_gemellarii"      "Lampanyctus_crocodilus"    
[15] "Holtbyrnia_macrops"         "Ceratoscopelus_maderensis" 
[17] "Chauliodus_sloani"          "Evermannella_balbo"        
[19] "Maulisia_mauli"             "Derichthys_serpentinus"    
[21] "Paralepis_coregonoides"     "Argyropelecus_hemigymnus"  
[23] "Melanostomias_bartonbeani" 

$Bathypelagic
 [1] "Argyropelecus_olfersii"    "Paralepis_coregonoides"   
 [3] "Borostomias_antarcticus"   "Malacosteus_niger"        
 [5] "Melanostigma_atlanticum"   "Argyropelecus_hemigymnus" 
 [7] "Gonostoma_elongatum"       "Arctozenus_risso"         
 [9] "Serrivomer_beanii"         "Maulisia_microlepis"      
[11] "Stomias_boa"               "Lestidiops_sphyrenoides"  
[13] "Bathylagus_euryops"        "Sigmops_bathyphilus"      
[15] "Lampanyctus_crocodilus"    "Maulisia_argipalla"       
[17] "Ceratoscopelus_maderensis" "Maulisia_mauli"           
[19] "Benthosema_glaciale"       "Evermannella_balbo"       
[21] "Chauliodus_sloani"         "Maurolicus_muelleri"      
[23] "Lampanyctus_macdonaldi"    "Anoplogaster_cornuta"     
[25] "Derichthys_serpentinus"    "Holtbyrnia_anomala"       
[27] "Melanostomias_bartonbeani"

$Epipelagic
 [1] "Stomias_boa"               "Melanostigma_atlanticum"  
 [3] "Argyropelecus_olfersii"    "Arctozenus_risso"         
 [5] "Borostomias_antarcticus"   "Ceratoscopelus_maderensis"
 [7] "Xenodermichthys_copei"     "Lestidiops_sphyrenoides"  
 [9] "Lobianchia_gemellarii"     "Myctophum_punctatum"      
[11] "Maurolicus_muelleri"       "Benthosema_glaciale"      
[13] "Chauliodus_sloani"         "Paralepis_coregonoides"   
[15] "Lampanyctus_crocodilus"    "Serrivomer_beanii"        
[17] "Argyropelecus_hemigymnus"  "Melanostomias_bartonbeani"
[19] "Sigmops_bathyphilus"      

$`Lower mesopelagic`
 [1] "Stomias_boa"                 "Melanostigma_atlanticum"    
 [3] "Argyropelecus_olfersii"      "Xenodermichthys_copei"      
 [5] "Serrivomer_beanii"           "Borostomias_antarcticus"    
 [7] "Lestidiops_sphyrenoides"     "Gonostoma_elongatum"        
 [9] "Lampanyctus_crocodilus"      "Arctozenus_risso"           
[11] "Bolinichthys_supralateralis" "Maulisia_mauli"             
[13] "Benthosema_glaciale"         "Ceratoscopelus_maderensis"  
[15] "Holtbyrnia_macrops"          "Maurolicus_muelleri"        
[17] "Sagamichthys_schnakenbecki"  "Chauliodus_sloani"          
[19] "Evermannella_balbo"          "Paralepis_coregonoides"     
[21] "Maulisia_argipalla"          "Derichthys_serpentinus"     
[23] "Argyropelecus_hemigymnus"    "Melanostomias_bartonbeani"  

• overlap between convex hulls shaping each of the two species assemblages.

Code

beta_plot_fish <- mFD::beta.multidim.plot(
  output_beta_fd_multidim = beta_fd_indices_fish,
  plot_asb_nm             = c("Epipelagic", "Upper mesopelagic"),
  beta_family             = c("Jaccard"),
  plot_sp_nm              = c("Maurolicus_muelleri", "Lampanyctus_crocodilus", "Argyropelecus_olfersii"),
  faxes                   = paste0("PC", 1:4),
  name_file               = NULL,
  faxes_nm                = NULL,
  range_faxes             = c(NA, NA),
  check_input             = TRUE)

beta_plot_fish$"patchwork"

Code

beta_plot_fish <- mFD::beta.multidim.plot(
  output_beta_fd_multidim = beta_fd_indices_fish,
  plot_asb_nm             = c("Upper mesopelagic", "Lower mesopelagic"),
  beta_family             = c("Jaccard"),
  plot_sp_nm              = c("Maulisia_mauli", "Evermannella_balbo", "Borostomias_antarcticus"),
  faxes                   = paste0("PC", 1:4),
  name_file               = NULL,
  faxes_nm                = NULL,
  range_faxes             = c(NA, NA),
  check_input             = TRUE)

beta_plot_fish$"patchwork"

Code

beta_plot_fish <- mFD::beta.multidim.plot(
  output_beta_fd_multidim = beta_fd_indices_fish,
  plot_asb_nm             = c("Lower mesopelagic", "Bathypelagic"),
  beta_family             = c("Jaccard"),
  plot_sp_nm              = c("Anoplogaster_cornuta", "Holtbyrnia_anomala", "Photostylus_pycnopterus"),
  faxes                   = paste0("PC", 1:4),
  name_file               = NULL,
  faxes_nm                = NULL,
  range_faxes             = c(NA, NA),
  check_input             = TRUE)

beta_plot_fish$"patchwork"

Code

## Compute the range of functional axes:
range_sp_coord  <- range(sp_faxes_coord_fish)

## Based on the range of species coordinates values, compute a nice range ...
## ... for functional axes:
range_faxes <- range_sp_coord +
  c(-1, 1) * (range_sp_coord[2] - range_sp_coord[1]) * 0.05


####### Create a list that will contains plots for each combination of axis:
plot_FRic <- list()

####### Compute all the combiantion we can get and the number of plots
axes_plot <- utils::combn(c("PC1", "PC2", "PC3", "PC4"), 2)
plot_nb   <- ncol(axes_plot)


######## Loop on all pairs of axes:
# for each combinaison of two axis:
for (k in (1:plot_nb)) {
  
  # get names of axes to plot:
  xy_k <- axes_plot[1:2, k]
  
  
  ####### Steps previously showed
  
  # a - Background:
  # get species coordinates along the two studied axes:
  sp_faxes_coord_xy <- sp_faxes_coord_fish[, xy_k]
  
  # Plot background with grey backrgound:
  plot_k <- mFD::background.plot(range_faxes = range_faxes, 
                                 faxes_nm = c(xy_k[1], xy_k[2]),
                                 color_bg = "grey95")
  
  
  # b - Global convex-hull:
  # Retrieve vertices coordinates along the two studied functional axes:
  vert <- mFD::vertices(sp_faxes_coord = sp_faxes_coord_xy,  
                        order_2D = FALSE, 
                        check_input = TRUE)
  
  plot_k <- mFD::pool.plot(ggplot_bg = plot_k,
                           sp_coord2D = sp_faxes_coord_xy,
                           vertices_nD = vert,
                           plot_pool = FALSE,
                           color_pool = NA,
                           fill_pool = NA,
                           alpha_ch =  0.8,
                           color_ch = "white",
                           fill_ch = "white",
                           shape_pool = NA,
                           size_pool = NA,
                           shape_vert = NA,
                           size_vert = NA,
                           color_vert = NA,
                           fill_vert = NA)
  
  
  # c - Assemblages convex-hulls and species:
  
  # Step 1: Species coordinates:
  # Bathypelagic:
  ## filter species from Bathypelagic:
  sp_filter_Bathypelagic <- mFD::sp.filter(asb_nm = c("Bathypelagic"),
                                           sp_faxes_coord = sp_faxes_coord_xy,
                                           asb_sp_w = depth_fish_biomass)
  ## get species coordinates (Bathypelagic):
  sp_faxes_coord_Bathypelagic <- sp_filter_Bathypelagic$`species coordinates`
  
  # Lower mesopelagic:
  ## filter species from Lower mesopelagic:
  sp_filter_Lower_mesopelagic <- mFD::sp.filter(asb_nm = c("Lower mesopelagic"),
                                                sp_faxes_coord = sp_faxes_coord_xy,
                                                asb_sp_w = depth_fish_biomass)
  ## get species coordinates (Lower mesopelagic):
  sp_faxes_coord_Lower_mesopelagic <- sp_filter_Lower_mesopelagic$`species coordinates`
  
  # Upper mesopelagic:
  ## filter species from Upper mesopelagic:
  sp_filter_Upper_mesopelagic <- mFD::sp.filter(asb_nm = c("Upper mesopelagic"),
                                                sp_faxes_coord = sp_faxes_coord_xy,
                                                asb_sp_w = depth_fish_biomass)
  ## get species coordinates (Upper mesopelagic):
  sp_faxes_coord_Upper_mesopelagic <- sp_filter_Upper_mesopelagic$`species coordinates`
  
  # Epipelagic:
  ## filter species from Epipelagic:
  sp_filter_Epipelagic <- mFD::sp.filter(asb_nm = c("Epipelagic"),
                                         sp_faxes_coord = sp_faxes_coord_xy,
                                         asb_sp_w = depth_fish_biomass)
  ## get species coordinates (Epipelagic):
  sp_faxes_coord_Epipelagic <- sp_filter_Epipelagic$`species coordinates`
  
  # Step 1 follow-up Vertices names:
  # Bathypelagic:
  vert_nm_Bathypelagic <- mFD::vertices(sp_faxes_coord = sp_faxes_coord_Bathypelagic,
                                        order_2D = TRUE, 
                                        check_input = TRUE)
  
  # Lower mesopelagic:
  vert_nm_Lower_mesopelagic <- mFD::vertices(sp_faxes_coord = sp_faxes_coord_Lower_mesopelagic,
                                             order_2D = TRUE, 
                                             check_input = TRUE)
  
  # Upper mesopelagic:
  vert_nm_Upper_mesopelagic <- mFD::vertices(sp_faxes_coord = sp_faxes_coord_Upper_mesopelagic,
                                             order_2D = TRUE, 
                                             check_input = TRUE)
  
  # Epipelagic:
  vert_nm_Epipelagic <- mFD::vertices(sp_faxes_coord = sp_faxes_coord_Epipelagic,
                                      order_2D = TRUE, 
                                      check_input = TRUE)
  
  # Step 2: plot convex-hulls and species of studied assemblages:
  plot_k <- mFD::fric.plot(ggplot_bg = plot_k,
                           asb_sp_coord2D = list("Bathypelagic" = sp_faxes_coord_Bathypelagic,
                                                 "Lower mesopelagic" = sp_faxes_coord_Lower_mesopelagic,
                                                 "Upper mesopelagic" = sp_faxes_coord_Upper_mesopelagic,
                                                 "Epipelagic"= sp_faxes_coord_Epipelagic),
                           asb_vertices_nD = list("Bathypelagic" = vert_nm_Bathypelagic,
                                                  "Lower mesopelagic" = vert_nm_Lower_mesopelagic,
                                                  "Upper mesopelagic" = vert_nm_Upper_mesopelagic,
                                                  "Epipelagic"= vert_nm_Epipelagic),
                           
                           plot_sp = T,
                           
                           color_ch =  c("Bathypelagic" = "#3C685A",
                                         "Lower mesopelagic" = "#6255B4",
                                         "Upper mesopelagic" = "#D62246",
                                         "Epipelagic" = "#FEA520"),
                           fill_ch = c("Bathypelagic" = "#3C685A",
                                       "Lower mesopelagic" = "#6255B4",
                                       "Upper mesopelagic" = "#D62246",
                                       "Epipelagic" = "#FEA520"),
                           alpha_ch = c("Bathypelagic" = 0.2,
                                        "Lower mesopelagic" = 0.2,
                                        "Upper mesopelagic" = 0.2,
                                        "Epipelagic" = 0.2),
                           
                           shape_sp = c("Bathypelagic" = 20,
                                        "Lower mesopelagic" = 20,
                                        "Upper mesopelagic" = 20,
                                        "Epipelagic" = 20),
                           size_sp = c("Bathypelagic" = 0.4,
                                       "Lower mesopelagic" = 0.4,
                                       "Upper mesopelagic" = 0.4,
                                       "Epipelagic" = 0.4),
                           color_sp = c("Bathypelagic" = "#3C685A",
                                        "Lower mesopelagic" = "#6255B4",
                                        "Upper mesopelagic" = "#D62246",
                                        "Epipelagic" = "#FEA520"),
                           fill_sp = c("Bathypelagic" = "white",
                                       "Lower mesopelagic" = "#6255B4",
                                       "Upper mesopelagic" = "#D62246",
                                       "Epipelagic" = "#FEA520"),
                           
                           shape_vert = c("Bathypelagic" = 17,
                                          "Lower mesopelagic" = 17,
                                          "Upper mesopelagic" = 17,
                                          "Epipelagic" = 17),
                           size_vert = c("Bathypelagic" = 4,
                                         "Lower mesopelagic" = 4,
                                         "Upper mesopelagic" = 4,
                                         "Epipelagic" = 4),
                           color_vert = c("Bathypelagic" = "#3C685A",
                                          "Lower mesopelagic" = "#6255B4",
                                          "Upper mesopelagic" = "#D62246",
                                          "Epipelagic" = "#FEA520"),
                           fill_vert = c("Bathypelagic" = "white",
                                         "Lower mesopelagic" = "#6255B4",
                                         "Upper mesopelagic" = "#D62246",
                                         "Epipelagic" = "#FEA520"))
  ####### Save the plot on the plot list:
  plot_FRic[[k]] <- plot_k
  
}

#plot_FRic

patchwork_FRic <- (plot_FRic[[1]] + patchwork::plot_spacer() + patchwork::plot_spacer() +
                     plot_FRic[[2]] + plot_FRic[[4]] + patchwork::plot_spacer() +
                     plot_FRic[[3]] + plot_FRic[[5]] + plot_FRic[[6]]) +
  patchwork::plot_layout(byrow = TRUE, heights = rep(1, 3),
                         widths = rep(1, 3), ncol = 3, nrow = 3,
                         guides = "collect")
patchwork_FRic

Code

ggsave("beta_diversity.png", path = "figures/additional_analyses", height = 9, width = 10)

\(\alpha\) & \(\beta\) diversity

• Moderate species turnover (29 %) and very low functional turnover (1% ).
• Most of the trait dissimilarity between species is found within a depth layer and not across depths

Code

#  Alpha, Beta and Gamma FD----
source(here::here("R", "Rao.r"))

partRao <-
  Rao(
    spxcom.matrix,
    dfunc = FD::gowdis(as.data.frame(spxtraits.matrix)),
    dphyl = NULL,
    weight = F,
    Jost = T,
    structure = NULL
  )

# Calculate proportions
result_df <- data.frame(
  Category = rep(c("Functional diversity", "Taxonomic diversity"), each = 2),
  Metric = c("Alpha", "Beta"),
  Value = c(
    100 * partRao$FD$Mean_Alpha / partRao$FD$Gamma,
    partRao$FD$Beta_prop,
    100 * partRao$TD$Mean_Alpha / partRao$TD$Gamma,
    partRao$TD$Beta_prop
  )
)

result_df$Category <- factor(result_df$Category, levels = c("Taxonomic diversity", "Functional diversity"))

ggplot(result_df, aes(x = Category, y = Value, fill = Metric)) +
  geom_bar(position = "stack", stat = "identity", alpha = 0.6) +
  scale_fill_manual(values = c("#338BA7", "#F59432")) +
  labs(x = "", y = "Proportion of Alpha and Beta in Equivalent numbers", fill = "") +
  geom_text(aes(label = sprintf("%.2f", Value)), vjust = -0.5, position = position_stack(vjust = 0.5)) +
  theme_minimal()

Code

round(partRao$TD$Pairwise_samples$Beta_prop / (1 - 1 / 3))

                  Epipelagic Upper mesopelagic Lower mesopelagic
Upper mesopelagic         41                                    
Lower mesopelagic         33                11                  
Bathypelagic              31                39                25

3.3. Functional redundancy

• Species Richness (N)
• Simpson diversity (D)
• Rao diversity (Q) : Sum of distances between pairs of randomly chosen species in trait space weighted by a relative abundance
• Functional redundancy : difference between species diversity and Rao’s quadratic entropy based on their functional dissimilarity (de Bello et al. 2007)
• Functionnal redundancy Ricotta = 1-Q/D
• Comparison with null models taking species diversity into account ? (SESRoa, SESRedon…)

Code

#Depth biomass data 
depth_fish_biomass_indices <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered, depth)%>%
  # divise biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
  select(species, depth, biomass_cpu)%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_cpu))%>%
  select(-c(biomass_cpu))%>%
  distinct()%>%
  arrange(depth) %>% 
  tidyr::pivot_wider(names_from = species, values_from = biomass)%>%
  replace(is.na(.), 0)%>%
  tibble::column_to_rownames(var = "depth")%>%
  #change species name
  rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
  rename("Cyclothone_sp"="Cyclothone")%>%
  rename("Stomias_boa"="Stomias_boa_boa") %>%
  as.matrix()

# Species richness ----
nsp <- as.data.frame(vegan::specnumber(depth_fish_biomass_indices)) %>%
  tibble::rownames_to_column(var = "depth") %>% 
  mutate(depth=as.numeric(depth))
colnames(nsp) <- c("depth", "species_richness")

# Diversity indices
diversity_indices_list <- SYNCSA::rao.diversity(depth_fish_biomass_indices, traits = fish_traits) 

diversity_indices <- data.frame(diversity_indices_list[c(2:4)]) %>% 
  tibble::rownames_to_column(var = "depth") %>% 
  mutate(depth=as.numeric(depth)) %>% 
  inner_join(nsp) %>% 
  mutate(functionnal_redundancy_ricotta= 1-(FunRao/Simpson)) %>% 
  tidyr::pivot_longer(! depth, names_to = "indices") %>% 
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

diversity_indices$depth_layer <-
  factor(
    diversity_indices$depth_layer,
    levels = c(
      "Epipelagic",
      "Upper mesopelagic",
      "Lower mesopelagic",
      "Bathypelagic"
    )
  )

diversity_indices$indices <- factor(
  diversity_indices$indices,
  levels = c(
    "species_richness",
    "Simpson",
    "FunRao",
    "FunRedundancy",
    "functionnal_redundancy_ricotta"
  ),
  labels = c(
    "Species richness (N)",
    "Simpson diversity (D)",
    "Rao diversity (Q)",
    "Functional Redundancy (R)",
    "Functional Redundancy (R=1-(Q/D))"
  )
)

ggplot(diversity_indices, aes(x = depth, y=value))+
  geom_point(alpha=0.5, size=1.5)+
  geom_smooth(method = "lm", col="#008E9B", alpha=0.1)+
  facet_wrap(~indices, scales = "free")+
  ggpmisc::stat_poly_eq(formula = y ~ x, 
                        aes(label = paste(..rr.label.., ..p.value.label.. 
                                          , ..n.label..,sep = "*`,`~")),
                        parse = TRUE,
                        size=3.7,
                        label.x.npc = "right",
                        label.y.npc = "bottom",
                        vstep = -0.0005)+
  theme_light()+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text(size =14, face="bold"),  
        legend.title = element_text(size =11),  
        legend.text = element_text(size =11), 
        axis.title.y = element_text(size=11),
        axis.text.y = element_text(size=12),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 11, face = "bold", color = "black"))

Code

ggsave("diversity_indices_lm.png", path = "figures/additional_analyses", height = 7, width = 10)

Code

ggplot(diversity_indices, aes(x =depth_layer, y=value))+
  geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.1) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  facet_wrap(~indices, scales = "free")+
  theme_light()+
  guides(col="none", fill="none")+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text(size =14, face="bold"),  
        legend.title = element_text(size =11),  
        legend.text = element_text(size =11), 
        axis.title.y = element_text(size=11),
        axis.text.y = element_text(size=12),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 11, face = "bold", color = "black"))

Code

ggsave("diversity_indices.png", path = "figures/additional_analyses", height = 7, width = 12)

4. Functional rarity

4.1 Different indices of functional rarity

Functional originality indices:

• Functional distinctiveness is the mean of dissimilarity of the focal species to all the other species of the set of interest. It can be abundance-weighted if needed.

• Functional uniqueness is the smallest dissimilarity that exists between the focal species and the all other species in the set. It does not consider the abundance of any species.

Rarity indices:

• Scarcity is proportional to the relative abundance of the species. It gets close to one when the species is (relatively) rare and close to 0 when its dominant

• Restrictedness is 1 minus the ratio of sites a species occupy over the total number of sites.

4.2.Computing functional rarity

4.2.1 Functional originality at regional scale

• For the choice or dissimilarity matrix we can use the raw dissimilarity matrix computed directly on raw traits values among species:

Code

sp_di <- funrar::distinctiveness_global(sp_dist_fish, di_name = "distinctiveness")

htmltools::tagList(DT::datatable(sp_di))

To compute uniqueness at regional scale we also need the regional level functional dissimilarity matrix with the uniqueness() function, and the site-species matrix:

Code

sp_ui <- funrar::uniqueness(
  pres_matrix = depth_fish_biomass,
  as.matrix(sp_dist_fish)
)

quantile(sp_ui$Ui, probs = seq(0, 1, by = 0.1))

        0%        10%        20%        30%        40%        50%        60% 
0.04008996 0.04907885 0.05443862 0.07016673 0.07846449 0.11224542 0.12990864 
       70%        80%        90%       100% 
0.14476124 0.17153014 0.19813792 0.23846021 

Code

htmltools::tagList(DT::datatable(sp_ui))

Based on these results we see that Anoplogaster cornuta, and Malacosteus niger are the most isolated fish in the functional space. Meaning that they have the most distant nearest neighbors.

4.2.2 Functional originality at local scale

Code

sp_local_di <- funrar::distinctiveness(
  depth_fish_biomass, as.matrix(sp_dist_fish)
)
sp_local_di[1:4, 1:10]

                  Anoplogaster_cornuta Arctozenus_risso
Upper mesopelagic                   NA        0.2270179
Bathypelagic                 0.3910062        0.2310470
Epipelagic                          NA        0.2170128
Lower mesopelagic                   NA        0.2281795
                  Argyropelecus_hemigymnus Argyropelecus_olfersii
Upper mesopelagic                0.2747268              0.2865240
Bathypelagic                     0.3025969              0.3054772
Epipelagic                       0.2907278              0.3170156
Lower mesopelagic                0.2919239              0.2906981
                  Bathylagus_euryops Benthosema_glaciale
Upper mesopelagic                 NA           0.2278397
Bathypelagic                0.215882           0.2107674
Epipelagic                        NA           0.1825419
Lower mesopelagic                 NA           0.2201022
                  Bolinichthys_supralateralis Borostomias_antarcticus
Upper mesopelagic                          NA               0.2907042
Bathypelagic                        0.2192975               0.2721636
Epipelagic                                 NA               0.2832503
Lower mesopelagic                   0.2383869               0.2791292
                  Ceratoscopelus_maderensis Chauliodus_sloani
Upper mesopelagic                 0.2069991         0.4064847
Bathypelagic                      0.1875186         0.4211923
Epipelagic                        0.1641562         0.4131897
Lower mesopelagic                 0.1976540         0.4039099

Code

identical(dim(sp_local_di), dim(depth_fish_biomass))

[1] TRUE

To compute uniqueness at the site scale, we must use a more complex expression as it was not envisioned for local computation:

Code

depth_ui <- apply(
  depth_fish_biomass, 1,
  function(single_site, dist_m) {
    single_site = single_site[single_site > 0 & !is.na(single_site)]
    funrar::uniqueness(t(as.matrix(single_site)), dist_m)
  }, dist_m = as.matrix(sp_dist_fish)
)

head(depth_ui[1])

$`Upper mesopelagic`
                      species         Ui
1            Arctozenus_risso 0.05200128
2    Argyropelecus_hemigymnus 0.11534810
3      Argyropelecus_olfersii 0.11534810
4         Benthosema_glaciale 0.07846449
5     Borostomias_antarcticus 0.14476124
6   Ceratoscopelus_maderensis 0.04907885
7           Chauliodus_sloani 0.22002144
8               Cyclothone_sp 0.20261783
9      Derichthys_serpentinus 0.17362662
10         Evermannella_balbo 0.23064259
11         Holtbyrnia_macrops 0.06075521
12     Lampanyctus_crocodilus 0.07016673
13    Lestidiops_sphyrenoides 0.05200128
14      Lobianchia_gemellarii 0.08729535
15          Malacosteus_niger 0.24719878
16             Maulisia_mauli 0.09014050
17        Maurolicus_muelleri 0.10554630
18    Melanostigma_atlanticum 0.14524258
19  Melanostomias_bartonbeani 0.19813792
20        Myctophum_punctatum 0.04907885
21           Lampanyctus_ater 0.07016673
22       Notoscopelus_kroyeri 0.05450362
23     Paralepis_coregonoides 0.14450490
24 Sagamichthys_schnakenbecki 0.06075521
25           Searsia_koefoedi 0.09014050
26          Serrivomer_beanii 0.17258393
27                Stomias_boa 0.19813792
28      Xenodermichthys_copei 0.14137149

As we had to manually build the function to compute the local uniqueness the results are strangely formatted.

We provide here a function that can help them to be more easily read:

Code

depth_ui <- lapply(names(depth_ui), function(x) {
  single_depth = depth_ui[[x]]
  single_depth$site = x
  
  return(single_depth)
})

depth_ui <- do.call(rbind, depth_ui)

#Then we can again look at the apple to see how its uniqueness varies across depths.

subset(depth_ui, species == "Melanostomias_bartonbeani")

                      species        Ui              site
19  Melanostomias_bartonbeani 0.1981379 Upper mesopelagic
56  Melanostomias_bartonbeani 0.1981379      Bathypelagic
83  Melanostomias_bartonbeani 0.1981379        Epipelagic
115 Melanostomias_bartonbeani 0.1981379 Lower mesopelagic

4.3.3.Rarity indices

scarcity:

Code

rel_weights = funrar::make_relative(depth_fish_biomass)

si =  funrar::scarcity(rel_weights)
summary(si)

 Anoplogaster_cornuta Arctozenus_risso Argyropelecus_hemigymnus
 Min.   :0.9307       Min.   :0.2023   Min.   :0.9565          
 1st Qu.:0.9307       1st Qu.:0.2305   1st Qu.:0.9806          
 Median :0.9307       Median :0.3366   Median :0.9901          
 Mean   :0.9307       Mean   :0.3458   Mean   :0.9826          
 3rd Qu.:0.9307       3rd Qu.:0.4520   3rd Qu.:0.9920          
 Max.   :0.9307       Max.   :0.5078   Max.   :0.9936          
 NA's   :3                                                     
 Argyropelecus_olfersii Bathylagus_euryops Benthosema_glaciale
 Min.   :0.1341         Min.   :0.1674     Min.   :0.09242    
 1st Qu.:0.4303         1st Qu.:0.1674     1st Qu.:0.10461    
 Median :0.5531         Median :0.1674     Median :0.34726    
 Mean   :0.4937         Mean   :0.1674     Mean   :0.42017    
 3rd Qu.:0.6166         3rd Qu.:0.1674     3rd Qu.:0.66282    
 Max.   :0.7345         Max.   :0.1674     Max.   :0.89375    
                        NA's   :3                             
 Bolinichthys_supralateralis Borostomias_antarcticus Ceratoscopelus_maderensis
 Min.   :0.9085              Min.   :0.7711          Min.   :0.1811           
 1st Qu.:0.9279              1st Qu.:0.9244          1st Qu.:0.2566           
 Median :0.9472              Median :0.9827          Median :0.4640           
 Mean   :0.9472              Mean   :0.9321          Mean   :0.4547           
 3rd Qu.:0.9666              3rd Qu.:0.9903          3rd Qu.:0.6622           
 Max.   :0.9859              Max.   :0.9919          Max.   :0.7098           
 NA's   :2                                                                    
 Chauliodus_sloani Cyclothone_sp    Derichthys_serpentinus
 Min.   :0.6243    Min.   :0.4460   Min.   :0.9488        
 1st Qu.:0.8282    1st Qu.:0.7333   1st Qu.:0.9633        
 Median :0.9415    Median :0.8690   Median :0.9779        
 Mean   :0.8768    Mean   :0.7822   Mean   :0.9710        
 3rd Qu.:0.9901    3rd Qu.:0.9179   3rd Qu.:0.9821        
 Max.   :1.0000    Max.   :0.9448   Max.   :0.9864        
                                    NA's   :1             
 Diaphus_metopoclampus Evermannella_balbo Gonostoma_elongatum
 Min.   :0.9889        Min.   :0.8112     Min.   :0.8448     
 1st Qu.:0.9904        1st Qu.:0.9025     1st Qu.:0.8754     
 Median :0.9920        Median :0.9939     Median :0.9060     
 Mean   :0.9920        Mean   :0.9340     Mean   :0.9060     
 3rd Qu.:0.9936        3rd Qu.:0.9954     3rd Qu.:0.9365     
 Max.   :0.9951        Max.   :0.9970     Max.   :0.9671     
 NA's   :2             NA's   :1          NA's   :2          
 Holtbyrnia_anomala Holtbyrnia_macrops Lampanyctus_crocodilus
 Min.   :0.8722     Min.   :0.9714     Min.   :0.003698      
 1st Qu.:0.8722     1st Qu.:0.9745     1st Qu.:0.006287      
 Median :0.8722     Median :0.9776     Median :0.022811      
 Mean   :0.8722     Mean   :0.9784     Mean   :0.048458      
 3rd Qu.:0.8722     3rd Qu.:0.9819     3rd Qu.:0.064982      
 Max.   :0.8722     Max.   :0.9862     Max.   :0.144510      
 NA's   :3          NA's   :1                                
 Lampanyctus_macdonaldi Lestidiops_sphyrenoides Lobianchia_gemellarii
 Min.   :0.3621         Min.   :0.9384          Min.   :0.7177       
 1st Qu.:0.3621         1st Qu.:0.9809          1st Qu.:0.8512       
 Median :0.3621         Median :0.9974          Median :0.9049       
 Mean   :0.3621         Mean   :0.9833          Mean   :0.8692       
 3rd Qu.:0.3621         3rd Qu.:0.9999          3rd Qu.:0.9229       
 Max.   :0.3621         Max.   :1.0000          Max.   :0.9493       
 NA's   :3                                                           
 Malacosteus_niger Maulisia_argipalla Maulisia_mauli   Maulisia_microlepis
 Min.   :0.8589    Min.   :0.9275     Min.   :0.4351   Min.   :0.5511     
 1st Qu.:0.8847    1st Qu.:0.9368     1st Qu.:0.7010   1st Qu.:0.5511     
 Median :0.9106    Median :0.9461     Median :0.9669   Median :0.5511     
 Mean   :0.9106    Mean   :0.9461     Mean   :0.7991   Mean   :0.5511     
 3rd Qu.:0.9364    3rd Qu.:0.9554     3rd Qu.:0.9812   3rd Qu.:0.5511     
 Max.   :0.9623    Max.   :0.9647     Max.   :0.9954   Max.   :0.5511     
 NA's   :2         NA's   :2          NA's   :1        NA's   :3          
 Maurolicus_muelleri Melanostigma_atlanticum Melanostomias_bartonbeani
 Min.   :0.1681      Min.   :0.9270          Min.   :0.9145           
 1st Qu.:0.7054      1st Qu.:0.9636          1st Qu.:0.9152           
 Median :0.9235      Median :0.9793          Median :0.9265           
 Mean   :0.7501      Mean   :0.9678          Mean   :0.9419           
 3rd Qu.:0.9681      3rd Qu.:0.9835          3rd Qu.:0.9532           
 Max.   :0.9853      Max.   :0.9855          Max.   :1.0000           
                                                                      
 Myctophum_punctatum Lampanyctus_ater Normichthys_operosus Notoscopelus_kroyeri
 Min.   :0.02587     Min.   :0.4081   Min.   :0.02368      Min.   :0.1825      
 1st Qu.:0.37286     1st Qu.:0.6404   1st Qu.:0.02368      1st Qu.:0.2149      
 Median :0.49108     Median :0.8437   Median :0.02368      Median :0.2863      
 Mean   :0.40536     Mean   :0.7678   Mean   :0.02368      Mean   :0.2769      
 3rd Qu.:0.52358     3rd Qu.:0.9711   3rd Qu.:0.02368      3rd Qu.:0.3483      
 Max.   :0.61339     Max.   :0.9758   Max.   :0.02368      Max.   :0.3523      
                                      NA's   :3                                
 Paralepis_coregonoides Photostylus_pycnopterus Sagamichthys_schnakenbecki
 Min.   :0.9793         Min.   :0.975           Min.   :0.9664            
 1st Qu.:0.9827         1st Qu.:0.975           1st Qu.:0.9721            
 Median :0.9859         Median :0.975           Median :0.9778            
 Mean   :0.9874         Mean   :0.975           Mean   :0.9751            
 3rd Qu.:0.9906         3rd Qu.:0.975           3rd Qu.:0.9795            
 Max.   :0.9983         Max.   :0.975           Max.   :0.9811            
                        NA's   :3               NA's   :1                 
 Searsia_koefoedi Serrivomer_beanii Sigmops_bathyphilus  Stomias_boa    
 Min.   :0.5166   Min.   :0.05462   Min.   :0.8561      Min.   :0.1152  
 1st Qu.:0.5425   1st Qu.:0.30778   1st Qu.:0.8880      1st Qu.:0.1431  
 Median :0.7013   Median :0.40401   Median :0.9199      Median :0.1813  
 Mean   :0.6953   Mean   :0.37437   Mean   :0.9199      Mean   :0.1968  
 3rd Qu.:0.8541   3rd Qu.:0.47060   3rd Qu.:0.9519      3rd Qu.:0.2350  
 Max.   :0.8618   Max.   :0.63484   Max.   :0.9838      Max.   :0.3093  
                                    NA's   :2                           
 Xenodermichthys_copei Notoscopelus_bolini
 Min.   :0.001126      Min.   :0.9986     
 1st Qu.:0.007770      1st Qu.:0.9993     
 Median :0.129474      Median :1.0000     
 Mean   :0.180620      Mean   :0.9995     
 3rd Qu.:0.302324      3rd Qu.:1.0000     
 Max.   :0.462404      Max.   :1.0000     
                       NA's   :1          

restrictiveness:

Code

ri = funrar::restrictedness(depth_fish_biomass)
summary(ri)

   species                Ri        
 Length:41          Min.   :0.0000  
 Class :character   1st Qu.:0.0000  
 Mode  :character   Median :0.0000  
                    Mean   :0.2378  
                    3rd Qu.:0.5000  
                    Max.   :0.7500  

4.3 Plotting functional rarity

4.3.1 Plotting functional originality

option to be able to colour species according to their functional originality (and not use the ready-made functions in the mfd package)

Code

# Make a summary data.frame
sp_coord_di_ui <- as.data.frame(sp_faxes_coord_fish[, 1:2])
sp_coord_di_ui$species <- rownames(sp_coord_di_ui)
rownames(sp_coord_di_ui) <- NULL
sp_coord_di_ui <- sp_coord_di_ui[, c(3, 1, 2)]
sp_coord_di_ui <- merge(sp_coord_di_ui, sp_di, by = "species")
sp_coord_di_ui <- merge(sp_coord_di_ui, sp_ui, by = "species")


plot_reg_distinctiveness <- ggplot(sp_coord_di_ui, aes(PC1, PC2)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point(aes(color = distinctiveness), size=3) +
  ggrepel::geom_text_repel(aes(label = species), size=4) +
  scale_color_gradient(high = "#914568", low = "#6BA1B9", "Functional\nDistinctiveness")+
  theme_bw() +
  theme(aspect.ratio = 1,
        legend.title = element_text(size=12),
        legend.text = element_text(size=12))

plot_reg_uniqueness <- ggplot(sp_coord_di_ui, aes(PC1, PC2)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point(aes(color = Ui), size=3) +
  ggrepel::geom_text_repel(aes(label = species), size=4) +
  scale_color_gradient(high = "#914568", low = "#6BA1B9","Functional\nUniqueness") +
  theme_bw() +
  theme(aspect.ratio = 1,
        legend.title = element_text(size=12),
        legend.text = element_text(size=12))

patchwork::wrap_plots(plot_reg_distinctiveness, plot_reg_uniqueness)

Code

ggsave("plot_reg_uniqueness.png", path = "figures", height = 12, width = 13)

Code

# Make a summary data.frame
sp_coord_di_ui <- as.data.frame(sp_faxes_coord_fish[, 1:2])
sp_coord_di_ui$species <- rownames(sp_coord_di_ui)
rownames(sp_coord_di_ui) <- NULL
sp_coord_di_ui <- sp_coord_di_ui[, c(3, 1, 2)]
sp_coord_di_ui <- merge(sp_coord_di_ui, sp_di, by = "species")
sp_coord_di_ui <- merge(sp_coord_di_ui, sp_ui, by = "species")
sp_coord_di_ui <- sp_coord_di_ui %>%
  merge(taxonomic_families, by = "species") %>%
  mutate(shape = case_when(
    family == "Myctophidae" ~ "A",
    family == "Platytroctidae" ~ "B",
    TRUE ~ "C"
  ))
  
plot_reg_uniqueness <- ggplot(sp_coord_di_ui, aes(PC1, PC2)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point(aes(color = Ui, shape = shape), size = 3) +
  ggrepel::geom_text_repel(aes(label = species), size = 3.5, max.overlaps = 5) +
  scale_color_viridis_c( name = "Functional Uniqueness", option="C") +
  theme_bw() +
  scale_shape_manual(values = c(15,17,19), name = "Taxonomic family",
                     labels = c("Myctophidae", "Platytroctidae", "Other"))+
  theme(aspect.ratio = 1, legend.title = element_text(size = 12), legend.text = element_text(size = 12))
plot_reg_uniqueness 

Code

ggsave("plot_reg_uniqueness.png", path = "figures", height = 8, width = 8)

• As was done with mFD to correlate the functional axes with species’ traits we can correlate functional distinctiveness to specific traits in order to see which traits are mainly driving distinctiveness.

• Regarding local level functional originality indices, the visualization can be more difficult to grasp and depends highly on the question. Would you rather focus on visualizing the functional distinctiveness of one species across communities? Compare the distribution of functional distinctiveness values across communities?

• One idea to keep in mind is that averaging functional distinctiveness per community is exactly equal to computing functional dispersion. Functional originality is computed on a species basis, so we should be aware that if we are rather interested by community properties than we can compute functional diversity metrics which are much more appropriate.

Code

local_di_ap <- as.data.frame(sp_local_di[, c(1: 11)])
local_di_ap$depth <- rownames(local_di_ap)

4.3.2. Plotting rarity

• Functional Distinctiveness

Code

depth_fish_biomass_2 <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered)%>%
  # divise biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
  select(species, depth, biomass_cpu)%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_cpu))%>%
  select(-c(biomass_cpu))%>%
  distinct()%>%
  tidyr::pivot_wider(names_from = species, values_from = biomass)%>%
  replace(is.na(.), 0)%>%
  tibble::column_to_rownames(var = "depth")%>%
  #change species name
  rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
  rename("Cyclothone_sp"="Cyclothone")%>%
  rename("Stomias_boa"="Stomias_boa_boa") %>%
  as.matrix()

ri = funrar::restrictedness(depth_fish_biomass_2)

median_depth_fish <- depth_distribution%>%
  mutate(species = case_when(
    species == "Nannobrachium_atrum"~"Lampanyctus_ater",
    species == "Cyclothone"~"Cyclothone_sp",
    species == "Stomias_boa_boa"~"Stomias_boa",
    TRUE ~ species
  )) %>% 
  group_by(species)%>%
  summarise(median_depth= median(depth))

sp_di_ri <- merge(sp_di, ri, by = "species")
sp_di_ri_ui <- merge(sp_di_ri, sp_ui, by = "species")%>%
  left_join(median_depth_fish)

ggplot(sp_di_ri, aes(distinctiveness, Ri)) +
  geom_point(aes(col=distinctiveness), size=3) +
  ggrepel::geom_text_repel(aes(label = species), max.overlaps=22, size=3) +
  labs(x = "Functional Distinctiveness", y = "Geographical Restrictedness") +
  theme_bw() +
  scale_color_viridis_c( name = "Functional Distinctiveness", option="C") +
  theme(aspect.ratio = 1,legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("sp_di_ri.png", path = "figures", height = 8, width = 8)

Code

ggplot(sp_di_ri_ui, aes(distinctiveness, Ri)) +
  geom_point(aes(col=median_depth), size=3) +
  ggrepel::geom_text_repel(aes(label = species), max.overlaps=22, size=3) +
  labs(x = "Functional Distinctiveness", y = "Geographical Restrictedness") +
  theme_bw() +
  scale_color_gradient(high = "#00263A", low = "#C5F8FF","Median depth (m)") +
  theme(aspect.ratio = 1,legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("sp_di_ri_depth.png", path = "figures", height = 8, width = 8)

• Functional uniqueness

Code

sp_di_ri <- merge(sp_di, ri, by = "species")
sp_di_ri_ui <- merge(sp_di_ri, sp_ui, by = "species")%>%
  left_join(median_depth_fish)

ggplot(sp_di_ri_ui, aes(Ui, Ri)) +
  geom_point(aes(col=median_depth), size=3) +
  ggrepel::geom_text_repel(aes(label = species), max.overlaps=22, size=3) +
  labs(x = "Functional Uniqueness", y = "Geographical Restrictedness") +
  theme_bw() +
  scale_color_gradient(high = "#00263A", low = "#C5F8FF","Median depth (m)") +
  theme(aspect.ratio = 1,legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("sp_ui_ri_depth.png", path = "figures", height = 8, width = 8)

• plot local scale measurements:

bathypelagic layer

Code

sp_local_di_df <- funrar::matrix_to_stack(
  sp_local_di, value_col = "local_di", row_to_col = "depth",
  col_to_col = "species"
)
sp_local_si_df <- funrar::matrix_to_stack(
  si, value_col = "local_si", row_to_col = "depth", col_to_col = "species"
)

sp_local_di_si <- merge(
  sp_local_di_df, sp_local_si_df, by = c("depth", "species")
)

head(sp_local_di_si)

         depth                  species  local_di  local_si
1 Bathypelagic     Anoplogaster_cornuta 0.3910062 0.9307396
2 Bathypelagic         Arctozenus_risso 0.2310470 0.5077535
3 Bathypelagic Argyropelecus_hemigymnus 0.3025969 0.9915088
4 Bathypelagic   Argyropelecus_olfersii 0.3054772 0.5289927
5 Bathypelagic       Bathylagus_euryops 0.2158820 0.1673770
6 Bathypelagic      Benthosema_glaciale 0.2107674 0.1086787

• At the local scale: Scarcity (relative abundance of the species) as a function of the local distinctiveness (one value per depth layer for each species)

Code

ggplot(sp_local_di_si, aes( local_si, local_di)) +
  geom_point(alpha = 0.4, size=1) +
  labs(x = "Functional Distinctiveness", y = "Scarcity") +
  theme_bw() +
  # geom_smooth(method = "lm",col="#9565E5", se=T, linewidth=1, fill="#9565E5", alpha=0.5)+
  # ggpmisc::stat_poly_eq(formula = y ~ x, 
  #                       aes(label = paste(..rr.label.., ..p.value.label.. 
  #                                         , ..n.label..,sep = "*`,`~")),
  #                       parse = TRUE,
  #                       size=2.5,
  #                       label.x.npc = "right",
  #                       label.y.npc = "bottom",
  #                       vstep = -0.0005)+
  # geom_text(aes(label = species), hjust = 0.7, vjust = -1, size = 3.5) +
  theme(aspect.ratio = 1)

Code

ggsave("local_di_si.png", path = "figures/additional_analyses", height = 5.5, width=5.5)

Code

ggplot(sp_di_ri_ui, aes(Ri,Ui)) +
  geom_point(alpha=0.6) +
  #ggrepel::geom_text_repel(aes(label = species), max.overlaps=22, size=3) +
  labs(x = "Functional Uniqueness", y = "Geographical Restrictedness") +
  theme_bw() +
  theme(aspect.ratio = 1)

Code

ggsave("global_ri_ui.png", path = "figures/additional_analyses", height = 5.5, width=5.5)

number families by depth layer

Code

 depth_fish_biomass %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column(var="assemblage") %>% 
  tidyr::pivot_longer(!assemblage, names_to = "species", values_to = "biomass") %>% 
  left_join(taxonomic_families) %>% 
  filter(biomass>0) %>% 
  select(assemblage, family) %>% 
   distinct() %>% 
  group_by(assemblage) %>% 
  summarize(n=n())

# A tibble: 4 × 2
  assemblage            n
  <chr>             <int>
1 Bathypelagic         14
2 Epipelagic           10
3 Lower mesopelagic    12
4 Upper mesopelagic    12

which traits are mainly driving distinctiveness?

Code

df <- fish_traits %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "species") %>%
  left_join(sp_di) %>%
  tibble::column_to_rownames(var = "species")

# Identify numerical and categorical traits
numeric_traits <-
  colnames(df)[sapply(df, is.numeric) &
                 colnames(df) != "distinctiveness"]
categorical_traits <-
  colnames(df)[!sapply(df, is.numeric) &
                 colnames(df) != "distinctiveness"]

# Initialize a data frame to store all results
combined_results_df <- data.frame(
  Trait = character(0),
  Type = character(0),
  Eta_R_squared = numeric(0),
  P_value = numeric(0),
  stringsAsFactors = FALSE
)

# Step 1: Perform Kruskal-Wallis test for categorical traits
for (categorical_trait in categorical_traits) {
  kruskal_result <-
    kruskal.test(df$distinctiveness ~ df[[categorical_trait]])
  
  # Calculate eta-squared for Kruskal-Wallis
  n_groups <- length(unique(df[[categorical_trait]]))
  n_total <- length(df$distinctiveness)
  h_value <- kruskal_result$statistic
  eta_squared <- (h_value - (n_groups - 1)) / (n_total - n_groups)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = categorical_trait,
      Type = "Categorical",
      Eta_R_squared = as.numeric(format(eta_squared, scientific = FALSE)),
      P_value = as.numeric(format(kruskal_result$p.value, scientific = FALSE)),
      stringsAsFactors = FALSE
    )
  )
}

# Step 2: Fit linear models for numerical traits
for (numeric_trait in numeric_traits) {
  lm_result <- lm(distinctiveness ~ df[[numeric_trait]], data = df)
  summary_stats <- summary(lm_result)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = numeric_trait,
      Type = "Numeric",
      Eta_R_squared = as.numeric(format(summary_stats$r.squared, scientific = FALSE)),
      P_value = as.numeric(format(
        summary_stats$coefficients[2, 4], scientific = FALSE
      )),
      stringsAsFactors = FALSE
    )
  )
}

# Round the numeric columns to three decimal places
combined_results_df[, c("Eta_R_squared", "P_value")] <-
  round(combined_results_df[, c("Eta_R_squared", "P_value")], 3)

htmltools::tagList(DT::datatable(combined_results_df))

which traits are mainly driving uniqueness ?

Code

df <- fish_traits %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "species") %>%
  left_join(sp_ui) %>%
  tibble::column_to_rownames(var = "species")

# Identify numerical and categorical traits
numeric_traits <-
  colnames(df)[sapply(df, is.numeric) &
                 colnames(df) != "Ui"]
categorical_traits <-
  colnames(df)[!sapply(df, is.numeric) &
                 colnames(df) != "Ui"]

# Initialize a data frame to store all results
combined_results_df <- data.frame(
  Trait = character(0),
  Type = character(0),
  Eta_R_squared = numeric(0),
  P_value = numeric(0),
  stringsAsFactors = FALSE
)

# Step 1: Perform Kruskal-Wallis test for categorical traits
for (categorical_trait in categorical_traits) {
  kruskal_result <-
    kruskal.test(df$Ui ~ df[[categorical_trait]])
  
  # Calculate eta-squared for Kruskal-Wallis
  n_groups <- length(unique(df[[categorical_trait]]))
  n_total <- length(df$Ui)
  h_value <- kruskal_result$statistic
  eta_squared <- (h_value - (n_groups - 1)) / (n_total - n_groups)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = categorical_trait,
      Type = "Categorical",
      Eta_R_squared = as.numeric(format(eta_squared, scientific = FALSE)),
      P_value = as.numeric(format(kruskal_result$p.value, scientific = FALSE)),
      stringsAsFactors = FALSE
    )
  )
}

# Step 2: Fit linear models for numerical traits
for (numeric_trait in numeric_traits) {
  lm_result <- lm(Ui ~ df[[numeric_trait]], data = df)
  summary_stats <- summary(lm_result)
  
  combined_results_df <- rbind(
    combined_results_df,
    data.frame(
      Trait = numeric_trait,
      Type = "Numeric",
      Eta_R_squared = as.numeric(format(summary_stats$r.squared, scientific = FALSE)),
      P_value = as.numeric(format(
        summary_stats$coefficients[2, 4], scientific = FALSE
      )),
      stringsAsFactors = FALSE
    )
  )
}

# Round the numeric columns to three decimal places
combined_results_df[, c("Eta_R_squared", "P_value")] <-
  round(combined_results_df[, c("Eta_R_squared", "P_value")], 3)

htmltools::tagList(DT::datatable(combined_results_df))

Code

df_long <- tidyr::gather(df, key = "Trait", value = "TraitValue", -Ui)

# Filter only significant traits
significant_traits <- combined_results_df[combined_results_df$P_value < 0.05, "Trait"]

df_long_filtered <- df_long[df_long$Trait %in% significant_traits, ]

# numeric traits 
numeric_traits <- df_long_filtered %>% 
  filter(Trait%in% c("head_length","eye_size", 
                     "pectoral_fin_insertion",
                     "eye_position", "orbital_length"))

numeric_traits$TraitValue <- as.numeric(numeric_traits$TraitValue)

ggplot(numeric_traits, aes(x = TraitValue, y = Ui)) +
  labs(x = "Trait Value", y = "Uniqueness") +
  theme_minimal() +
  facet_wrap(~ Trait, scales = "free") +
  geom_point() +
  geom_smooth(method = "lm") 

Code

categorical_traits <- df_long_filtered %>% 
  filter(!Trait%in% c("head_length","eye_size", 
                     "pectoral_fin_insertion",
                     "eye_position", "orbital_length"))

ggplot(categorical_traits, aes(x = TraitValue, y = Ui)) +
  labs(x = "Trait Value", y = "Uniqueness") +
  theme_minimal() +
  facet_wrap(~ Trait, scales = "free") +
  geom_point() +
  geom_boxplot()

Relation between biomass and distinctiveness

Code

depth_biomass <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered)%>%
  # divise biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
  select(species, depth, biomass_cpu)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_cpu))%>%
  select(-c(biomass_cpu))%>%
  distinct()  %>%
  mutate(species = case_when(
    species == "Nannobrachium_atrum" ~ "Lampanyctus_ater",
    species == "Cyclothone" ~ "Cyclothone_sp",
    species == "Stomias_boa_boa" ~ "Stomias_boa",
    TRUE ~ species
  ))

Uniqueness - biomass log

Code

Ui_biomass <- depth_biomass %>% 
  group_by(species) %>% 
  replace(is.na(.), 0) %>% 
  mutate(biomass_sp = sum(biomass)) %>% 
  select(species, biomass_sp) %>% 
  distinct() %>% 
  inner_join(sp_ui)

ggplot(Ui_biomass , aes(x = Ui, y = log(biomass_sp))) +
  geom_point(size = 2, aes(col=Ui)) +
  scale_color_viridis_c( name = "Functional Uniqueness", option="C") +
  labs(x="Uniqueness", y= "Log biomass")+
  geom_text(aes(label = species), hjust = 0.7, vjust = -1, size = 3.5) +
  theme_bw()+
  theme(aspect.ratio = 1, legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("Ui_biomass_log.png", path = "figures", height = 8, width=8)

Uniqueness - abundance

Code

# list of species 
sp_names <- c(rownames(fish_traits), "Nannobrachium_atrum", "Cyclothone", "Stomias_boa_boa")

# Metadata
metadata <-  utils::read.csv(here::here("data", "metadata.csv"), sep = ";", header = T, dec = ".")%>%
  # calculation of standardized abundance values (vertical  trawl opening * horizontal trawl opening * distance traveled)  
  mutate(volume_filtered = 24*58*distance)

# species abundance x depth  matrix 2002-2019 ----
data_abundance_2002_2019 <- utils::read.csv(here::here("data", "data_evhoe_catch_2002_2019_abundance.csv"), sep = ";", header = T, dec = ".")%>%
  replace(is.na(.), 0)%>%
  as.data.frame()%>%
  rename("species"="Code_Station")%>%
  mutate(species= gsub(" ","_", species))%>%
  filter(species%in%sp_names)%>%
  t()%>%
  as.data.frame()%>%
  janitor::row_to_names(row_number = 1)%>%
  mutate_if(is.character, as.numeric)%>%
  tibble::rownames_to_column("Code_Station")%>%
  filter(!Code_Station=="H0472")%>%
  tidyr::pivot_longer(!Code_Station, names_to = "species", values_to = "Tot_V_HV")%>%
  rename("Nom_Scientifique"="species") %>% 
  rename("abundance"="Tot_V_HV")

# species abundance x depth  matrix 2021 ----
data_abundance_2021 <- utils::read.csv(here::here("data", "data_evhoe_catch_2021.csv"), sep = ";", header = T, dec = ".")%>%
  select(Nom_Scientifique, Nbr, Code_Station)%>%
  filter(Nbr>0) %>% 
  group_by(Nom_Scientifique, Code_Station) %>% 
  mutate(abundance=sum(Nbr)) %>% 
  select(Nom_Scientifique, abundance, Code_Station) %>% 
  distinct() %>% 
  mutate(Nom_Scientifique= gsub(" ","_",Nom_Scientifique))%>%
  filter(Nom_Scientifique%in%sp_names)

# species abundance x depth  matrix 2022 ----
data_abundance_2022 <- utils::read.csv(here::here("data", "data_evhoe_catch_2022.csv"), sep = ";", header = T, dec = ".")%>%
  select(Nom_Scientifique, Nbr, Code_Station)%>%
  filter(Nbr>0) %>% 
  group_by(Nom_Scientifique, Code_Station) %>% 
  mutate(abundance=sum(Nbr)) %>% 
  select(Nom_Scientifique, abundance, Code_Station) %>% 
  distinct() %>% 
  mutate(Nom_Scientifique= gsub(" ","_",Nom_Scientifique))%>%
  filter(Nom_Scientifique%in%sp_names)

# merge all matrix ----
depth_fish_abundance <- rbind(data_abundance_2002_2019, data_abundance_2021, data_abundance_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%  
  left_join(metadata) %>% 
  select(species, abundance, depth, volume_filtered)%>%
  # divise abundance by the volume filtered at each trawl (g.m3)
  mutate(abundance_cpu=(abundance/volume_filtered)*1000)%>%
  select(species, depth, abundance_cpu)%>%  
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  select(-depth)%>%
  group_by(species, depth_layer)%>%
  mutate(abundance=sum(abundance_cpu))%>%
  select(-c(abundance_cpu))%>%
  distinct()%>%
  tidyr::pivot_wider(names_from = species, values_from = abundance)%>%
  replace(is.na(.), 0)%>%
  tibble::column_to_rownames(var = "depth_layer")%>% 
  #change species name
  rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
  rename("Cyclothone_sp"="Cyclothone")%>%
  rename("Stomias_boa"="Stomias_boa_boa") %>% 
  as.matrix()

depth_abundance <- depth_fish_abundance %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "assemblage") %>%
  tidyr::pivot_longer(!assemblage, names_to = "species", values_to = "abundance") %>%
  filter(abundance > 0)

Ui_abundance <- depth_abundance %>% 
  group_by(species) %>% 
  mutate(abundance_sp = sum(abundance)) %>% 
  select(species, abundance_sp) %>% 
  distinct() %>% 
  inner_join(sp_ui)

ggplot(Ui_abundance , aes(x = Ui, y = log(abundance_sp))) +
  geom_point(size = 2, aes(col=Ui)) +
  scale_color_viridis_c( name = "Functional Uniqueness", option="C") +
  labs(x="Uniqueness")+
  geom_text(aes(label = species), hjust = 1, vjust = -1, size = 3.5) +
  theme_bw()+
  theme(aspect.ratio = 1, legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("Ui_abundance_log.png", path = "figures", height = 8, width=8)

Distinctiveness vs abundance

Code

Di_abundance <- depth_abundance %>% 
  group_by(species) %>% 
  mutate(abundance_sp = sum(abundance)) %>% 
  select(species, abundance_sp) %>% 
  distinct() %>% 
  inner_join(sp_di)

ggplot(Di_abundance , aes(x = distinctiveness, y = log(abundance_sp))) +
  geom_point(size = 2, aes(col=distinctiveness)) +
  scale_color_viridis_c( name = "Functional distinctiveness", option="C") +
  labs(x="Distinctiveness")+
  geom_text(aes(label = species), hjust = 1, vjust = -1, size = 3.5) +
  theme_bw()+
  theme(aspect.ratio = 1, legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("Di_abundance_log.png", path = "figures", height = 8, width=8)

Code

Di_biomass <- depth_biomass %>% 
  group_by(species) %>% 
  mutate(biomass_sp = sum(biomass)) %>% 
  select(species, biomass_sp) %>% 
  distinct() %>% 
  inner_join(sp_di)

ggplot(Di_biomass , aes(x = distinctiveness, y = log(biomass_sp))) +
  geom_point(size = 2, aes(col=distinctiveness)) +
  scale_color_viridis_c( name = "Functional distinctiveness", option="C") +
  labs(x="Distinctiveness")+
  geom_text(aes(label = species), hjust = 1, vjust = -1, size = 3.5) +
  theme_bw()+
  theme(aspect.ratio = 1, legend.title = element_text(size = 12), legend.text = element_text(size = 12))

Code

ggsave("Di_biomass_log.png", path = "figures", height = 8, width=8)